logo

Goal: append State Park footprints to preliminary CA-BBA blocks, visualize, and calculate basic summary statistics.

library(dplyr)
library(here)
library(mapview)
library(tigris)
library(sf)

mapviewOptions(basemaps.color.shuffle = FALSE)

ca_blocks_clean <- readRDS(here("Analysis", "Data", "Blocks", "v202506", "CA_blocks_v202506.RDS")) %>%
  filter(WATER_SLIVER == 0 & CA_SLIVER == 0)
cntr_crds <- c(mean(st_coordinates(ca_blocks_clean)[, 1]),
               mean(st_coordinates(ca_blocks_clean)[, 2]))

California State Parks boundaries as provided by Ron Melcer.

Multiple State Parks sometimes overlap a single block. In these cases, we record the ID of all parks but also determine (1) the “most-localized” park (the park with the greatest percentage of its total area located in the block) and (2) the “primary” park (the park with the largest footprint in the block). The map below shows each block’s most-localized park, but clicking individual blocks will also list the block’s primary park. Detailed park boundaries can be toggled on using the map’s layer button.

# Load State Park data
state_parks <- st_read(here("Analysis", "Data", "State_parks", "Park_Boundaries", "ParkBoundaries.shp")) %>%
  st_transform(st_crs(ca_blocks_clean)) %>%
  st_make_valid() %>%
  st_cast("POLYGON") %>%
  mutate("temp_ID" = row_number(.))

state_parks_temp <- state_parks %>%
  mutate(area = st_area(.) %>% units::drop_units()) %>%
  st_drop_geometry() %>%
  group_by(UNITNAME) %>%
  summarize(unitname_area =sum(area))

# Intersect with CA-BBA blocks, count number of parks by block, and record thier IDs and the ID of park with largest footprint in block
xFeatures <- ca_blocks_clean %>% select(BLOCKNAME)
yFeatures <- state_parks %>% select(temp_ID, FID, UNITNAME)
intersections <- st_intersects(x = xFeatures, y = yFeatures)

library(progress)
pb <- progress_bar$new(format = "[:bar] :current/:total (:percent)", total = dim(xFeatures)[1])

library(purrr)
state_park_intersect <- map_dfr(1:dim(xFeatures)[1], function(ix){
  pb$tick()
  st_intersection(x = xFeatures[ix,], y = yFeatures[intersections[[ix]],])
})

state_park_intersect_details <- state_park_intersect %>%
  mutate(area = st_area(geometry)) %>%
  as.data.frame() %>%
  st_drop_geometry() %>%
  select(BLOCKNAME, UNITNAME, area) %>%
  mutate(area = units::drop_units(area)) %>%
  group_by(BLOCKNAME, UNITNAME) %>%
  left_join(state_parks_temp %>% st_drop_geometry()) %>%
  summarize(park_area = sum(area),
            park_pct_area_in_block = (sum(area) / mean(unitname_area) * 100) %>% round(5),
            .groups = 'drop')
state_park_intersect <- state_park_intersect_details %>%
  group_by(BLOCKNAME) %>%
  summarize(
   STATE_PARK_COUNT = length(unique(UNITNAME)),
   STATE_PARK_PRIMARY = UNITNAME[which.max(park_area)],
   STATE_PARK_MOST_LOCALIZED = UNITNAME[which.max(park_pct_area_in_block)],
   STATE_PARK_MOST_LOCALIZED_PCT = park_pct_area_in_block[which.max(park_pct_area_in_block)],

   STATE_PARK_ALL = paste(unique(UNITNAME), collapse = "; "),
    .groups = 'drop'
  )

# Join State Park info back to full blocks
ca_blocks_clean <- ca_blocks_clean %>%
  left_join(state_park_intersect, by = "BLOCKNAME") 

# Calculate number of blocks with State Parks
n_blocks <- ca_blocks_clean %>% nrow()
n_blocks_SPs <- ca_blocks_clean %>% filter(!is.na(STATE_PARK_MOST_LOCALIZED)) %>% nrow()

Blocks in State Parks: 1050 of 16322 (6.4%)

 



By California Bird Atlas